This markdown provides comparison between two meQTL analysis after LD clumping. meQTL analysis were done using the R package MatrixEQTL. LD clumping was done using PLINK
For all analysis, the same DNAm data were used which underwent following QC steps:
DNAm QC
The genotype data were different for each analysis:
The SNP data were used where the uncertain genotype was replaced with the most likely genotype
SNP probabilities, i.e. dosages, were used
For all three analysis, the imputed QC SNP data (no LD pruning, the MHC regions are included) were used. The data includes 3’957’338 SNPs.
All meQTLs passed the significance threshold of FDR = 0.05.
LD clumping was done per CpG using the following parameters:
–clump-p1 = 0.05: Significance threshold for index SNPs
–clump-p2 = 1: Secondary significance threshold for clumped SNPs
–clump-r2 = 0.2: LD threshold for clumping
–clump-kb = 200: Physical distance threshold for clumping
x %>%
ggplot(aes(y = cnt, x = Var2, fill = as.factor(as.numeric(Var2)))) + geom_col() + geom_text(aes(label = comma(cnt,
accuracy = 1L), y = cnt), stat = "identity", vjust = -0.2, size = 2.5) + scale_y_continuous(labels = scientific) +
# scale_fill_viridis(discrete = TRUE, alpha = 0.9, option = 'C') +
labs(title = " ", y = "Count", x = " ") + facet_wrap(~Var1, ncol = 3) + theme(legend.position = "none",
legend.title = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(),
plot.title = element_text(size = 10), axis.title = element_text(size = 8), axis.title.x = element_blank()) +
scale_fill_manual(values = cbPalette[c(1, 2)])
fdr.thr <- 0.05
plot.title <- paste0("Significant at FDR <= ", fdr.thr, " meQTLs.", "\nResult of MatrixEQTL")
GetScatterPlot2(meqtl.all.prob.df, plot.title = plot.title, cbPalette = cbPalette)
meqtls.meqtls <- list(delta = meqtl.all.prob.df[treatment == "delta", meQTL_ID], veh = meqtl.all.prob.df[treatment ==
"baseline", meQTL_ID])
intersect.meqtls <- intersect(meqtls.meqtls$delta, meqtls.meqtls$veh)
perc.olap.delta.with.veh <- scales::percent(length(intersect.meqtls)/length(meqtls.meqtls$delta), accuracy = 0.1)
perc.olap.veh.with.delta <- scales::percent(length(intersect.meqtls)/length(meqtls.meqtls$veh), accuracy = 0.1)
plot.title <- "Number of intersections between delta and baseline meQTLs"
ggVennDiagram::ggVennDiagram(meqtls.meqtls, category.names = c(paste0("GR-response, ", perc.olap.delta.with.veh),
paste0("Baseline, ", perc.olap.veh.with.delta)), label_alpha = 0.7, edge_size = 0, set_geom = "text",
set_color = "black", label = "count") + theme(legend.position = "none", legend.title = element_blank(),
panel.background = element_blank(), plot.title = element_text(size = 10, face = "italic")) + labs(x = "",
y = "", title = plot.title) + scale_color_manual(values = cbPalette) + scale_fill_gradient(low = alpha(cbPalette[2],
0.5), high = alpha(cbPalette[1], 0.5))
meqtls.snps <- list(delta = meqtl.all.prob.df[treatment == "delta", SNP], veh = meqtl.all.prob.df[treatment ==
"baseline", SNP])
intersect <- intersect(meqtls.snps$delta, meqtls.snps$veh)
perc.olap.delta.with.veh <- scales::percent(length(intersect)/length(meqtls.snps$delta), accuracy = 0.1)
perc.olap.veh.with.delta <- scales::percent(length(intersect)/length(meqtls.snps$veh), accuracy = 0.1)
plot.title <- "Number of intersections between delta and baseline meQTL SNPs"
ggVennDiagram::ggVennDiagram(meqtls.snps, category.names = c(paste0("GR-response, ", perc.olap.delta.with.veh),
paste0("Baseline, ", perc.olap.veh.with.delta)), label_alpha = 0.7, edge_size = 0, set_geom = "text",
set_color = "black", label = "count") + theme(legend.position = "none", legend.title = element_blank(),
panel.background = element_blank(), plot.title = element_text(size = 10, face = "italic")) + labs(x = "",
y = "", title = plot.title) + scale_color_manual(values = cbPalette) + scale_fill_gradient(low = alpha(cbPalette[2],
0.5), high = alpha(cbPalette[1], 0.5))
meqtls.cpg <- list(delta = meqtl.all.prob.df[treatment == "delta", CpG_ID], veh = meqtl.all.prob.df[treatment ==
"baseline", CpG_ID])
intersect <- intersect(meqtls.cpg$delta, meqtls.cpg$veh)
perc.olap.delta.with.veh <- scales::percent(length(intersect)/length(meqtls.cpg$delta), accuracy = 0.1)
perc.olap.veh.with.delta <- scales::percent(length(intersect)/length(meqtls.cpg$veh), accuracy = 0.1)
plot.title <- "Number of intersections between delta and baseline meQTL CpGs"
ggVennDiagram::ggVennDiagram(meqtls.cpg, category.names = c(paste0("GR-response, ", perc.olap.delta.with.veh),
paste0("Baseline, ", perc.olap.veh.with.delta)), label_alpha = 0.7, edge_size = 0, set_geom = "text",
set_color = "black", label = "count") + theme(legend.position = "none", legend.title = element_blank(),
panel.background = element_blank(), plot.title = element_text(size = 10, face = "italic")) + labs(x = "",
y = "", title = plot.title) + scale_color_manual(values = cbPalette) + scale_fill_gradient(low = alpha(cbPalette[2],
0.5), high = alpha(cbPalette[1], 0.5))
x %>%
ggplot(aes(y = cnt, x = Var2, fill = as.factor(as.numeric(Var2)))) + geom_col() + geom_text(aes(label = comma(cnt,
accuracy = 1L), y = cnt), stat = "identity", vjust = -0.2, size = 2.5) + scale_y_continuous(labels = scientific) +
# scale_fill_viridis(discrete = TRUE, alpha = 0.9, option = 'C') +
labs(title = " ", y = "Count", x = " ") + facet_wrap(~Var1, ncol = 3) + theme(legend.position = "none",
legend.title = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(),
plot.title = element_text(size = 10), axis.title = element_text(size = 8), axis.title.x = element_blank()) +
scale_fill_manual(values = cbPalette[c(1, 2)])
meqtls.delta <- list(prob = meqtl.all.prob.df[treatment == "delta", meQTL_ID], mode = meqtl.all.mode.df[treatment ==
"delta", meQTL_ID])
meqtls.veh <- list(prob = meqtl.all.prob.df[treatment == "baseline", meQTL_ID], mode = meqtl.all.mode.df[treatment ==
"baseline", meQTL_ID])
euler.tbl <- euler(meqtls.delta)
plot(euler.tbl, quantities = list(type = c("counts", "percent")), fills = c("lightgrey", "brown"), labels = c("probabilities",
" most likely genotype"))
euler.tbl <- euler(meqtls.veh)
plot(euler.tbl, quantities = list(type = c("counts", "percent")), fills = c("lightgrey", "brown"), labels = c("probabilities",
" most likely genotype"))
mecpgs.delta <- list(prob = meqtl.all.prob.df[treatment == "delta", CpG_ID] %>%
unique(), mode = meqtl.all.mode.df[treatment == "delta", CpG_ID] %>%
unique())
mecpgs.veh <- list(prob = meqtl.all.prob.df[treatment == "baseline", CpG_ID] %>%
unique(), mode = meqtl.all.mode.df[treatment == "baseline", CpG_ID] %>%
unique())
euler.tbl <- euler(mecpgs.delta)
plot(euler.tbl, quantities = list(type = c("counts", "percent")), fills = c("lightgrey", "brown"), labels = c("probabilities",
" most likely genotype"))
euler.tbl <- euler(mecpgs.veh)
plot(euler.tbl, quantities = list(type = c("counts", "percent")), fills = c("lightgrey", "brown"), labels = c("probabilities",
" most likely genotype"))
mesnps.delta <- list(prob = meqtl.all.prob.df[treatment == "delta", SNP] %>%
unique(), mode = meqtl.all.mode.df[treatment == "delta", SNP] %>%
unique())
mesnps.veh <- list(prob = meqtl.all.prob.df[treatment == "baseline", SNP] %>%
unique(), mode = meqtl.all.mode.df[treatment == "baseline", SNP] %>%
unique())
euler.tbl <- euler(mesnps.delta)
plot(euler.tbl, quantities = list(type = c("counts", "percent")), fills = c("lightgrey", "brown"), labels = c("probabilities",
" most likely genotype"))
euler.tbl <- euler(mesnps.veh)
plot(euler.tbl, quantities = list(type = c("counts", "percent")), fills = c("lightgrey", "brown"), labels = c("probabilities",
" most likely genotype"))
meqtls.delta <- list(prob = meqtl.all.prob.df[treatment == "delta", meQTL_ID],
miss = NULL, # meqtl.all.na.df[treatment == "delta", meQTL_ID],
mode = meqtl.all.mode.df[treatment == "delta", meQTL_ID])
unique.mode.meqtl <- setdiff(meqtls.delta$mode, union(meqtls.delta$prob, meqtls.delta$miss))
unique.mode.meqtl.df <- meqtl.all.mode.df[treatment == "delta"][meQTL_ID %in% unique.mode.meqtl]
selected.meqtl <- unique.mode.meqtl.df[1]
kable(selected.meqtl %>%
select(SNP, CpG_ID, FC = beta, `t-stat`, `p-value`, FDR = fdr, Treatment = treatment) %>%
mutate_if(is.numeric, funs(as.character(signif(., 3)))))
| SNP | CpG_ID | FC | t-stat | p-value | FDR | Treatment |
|---|---|---|---|---|---|---|
| rs10249476 | cg00011113 | 0.0208 | 5.75 | 3.48e-08 | 8.72e-08 | delta |
ProcessGetBoxPlot(methyl.beta.veh.df, methyl.beta.dex.df, snp.df, selected.meqtl, fdr.thr = 0.05)
unique.meqtl <- setdiff(meqtls.delta$prob, union(meqtls.delta$mode, meqtls.delta$miss))
unique.meqtl.df <- meqtl.all.prob.df[treatment == "delta"][meQTL_ID %in% unique.meqtl]
selected.meqtl <- unique.meqtl.df[1]
kable(selected.meqtl %>%
select(SNP, CpG_ID, FC = beta, `t-stat`, `p-value`, FDR = fdr, Treatment = treatment) %>%
mutate_if(is.numeric, funs(as.character(signif(., 3)))))
| SNP | CpG_ID | FC | t-stat | p-value | FDR | Treatment |
|---|---|---|---|---|---|---|
| rs35879900 | cg26175789 | -0.031 | -5.44 | 1.62e-07 | 0.00875 | delta |
ProcessGetBoxPlot(methyl.beta.veh.df, methyl.beta.dex.df, snp.df, selected.meqtl, fdr.thr = 0.05)
meqtl.lst <- intersect(meqtls.delta$mode, meqtls.delta$prob)
meqtl.df <- meqtl.all.mode.df[treatment == "delta"][meQTL_ID %in% meqtl.lst]
selected.meqtl <- meqtl.df[1]
kable(selected.meqtl %>%
select(SNP, CpG_ID, FC = beta, `t-stat`, `p-value`, FDR = fdr, Treatment = treatment) %>%
mutate_if(is.numeric, funs(as.character(signif(., 3)))))
| SNP | CpG_ID | FC | t-stat | p-value | FDR | Treatment |
|---|---|---|---|---|---|---|
| rs114232786 | cg00020172 | 0.0452 | 6.41 | 1.06e-09 | 4.03e-09 | delta |
selected.meqtl <- meqtl.df[1]
ProcessGetBoxPlot(methyl.beta.veh.df, methyl.beta.dex.df, snp.df, selected.meqtl, fdr.thr = 0.05)
meqtl.df <- meqtl.all.prob.df[treatment == "delta"][meQTL_ID %in% meqtl.lst]
selected.meqtl <- meqtl.df[1]
kable(selected.meqtl %>%
select(SNP, CpG_ID, FC = beta, `t-stat`, `p-value`, FDR = fdr, Treatment = treatment) %>%
mutate_if(is.numeric, funs(as.character(signif(., 3)))))
| SNP | CpG_ID | FC | t-stat | p-value | FDR | Treatment |
|---|---|---|---|---|---|---|
| rs11107749 | cg08922308 | 0.0461 | 8 | 1.1e-13 | 4.93e-08 | delta |
selected.meqtl <- meqtl.df[1]
ProcessGetBoxPlot(methyl.beta.veh.df, methyl.beta.dex.df, snp.df, selected.meqtl, fdr.thr = 0.05)